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Features of the V3D-Vascular System 



Vessel displayed in curved reformat: A curved planar reformation view shall be 
available. 

o It shall take as input a 3D centerline. 

o It shall be possible to reformat the data in any of the following three ways: 

■ The most traditional method is what doctors routinely call the 
"curved MPR" or the "curved reformat". As shown above, the 
vessel will appear curved on the screen. If the vessel occludes 
itself (e.g., it loops around), only the front-most vessel section will 
be visible (i.e., it uses a z-buffer). 

■ The "luminal view" in which the vessel appears straight is an 
improved technique. In this technique, the scanlines of the image 
are acquired along the intersection of the image plane and the 
plane locally perpendicular to the centerline. 

o It shall be possible to rotate the angle of projection about the centerline 
360 degrees. 

■ The outline of the image forms a ribbon in 2D space and it should 
be displayed in correspondence on the 3D images. 

Vessel displayed in oblique reformat: An oblique planar reformation view shall 
be available. 

o The field of view shall be adjustable to accommodate the small coronary 

arteries up to the aorta, 
o It shall be possible to scroll along the centerline with the mouse wheel. 
Endoluminal view: An endo luminal view shall be available which makes the 
contrast material in the vessel (if any exists) transparent. 

o General interactive navigation through the vessel shall be possible at high 

frame rates (better than 2 frames per second), 
o The visualization shall enable the doctor to differentiate between different 

densities using color or some other scheme. This might enable the doctor 

to visualize soft plaque if it were sufficiently different in intensity from 

hard plaque and the average lumen density, 
o Guided navigation through the vessel from end to end shall be possible 

along the centerline. 
o An outline of the slice or slab shall be visible in correspondence on the 3D 

images. 

o Guided endoluminal flight shall be available from end to end and back, 
similar to the colon module. 
Volume rendering or MIP projection: The vascular dataset shall be viewable 
from the outside. 

o The projection shall either be using 3D volume rendering, X-ray, or MIP. 
o The view shall be rotatable and zoomable. 

o The image shall be labeled so that the six orthogonal directions are readily 
identifiable. 

Standard angiogram projections: There shall be a way to set up the projection 
along a variety of preset directions. 



o Once the projection direction is set up, it shall be possible to visualize 
using 3D rendering. 

Simplified selection of vessel trees: There shall be a specific segmentation mode 
that makes it easy (just a few clicks) to segment entire vessel trees if sufficient 
contrast is available. One click will be used to place a seedpoint, then a slider will 
be activated on the screen that only requires mouse movement (e.g., up/down) to 
control the amount of growth, then a second click to lock in the desired amount of 
growth. When selecting vessels, it will attempt to grow first through objects that 
look like vessels. When selecting non-vessels, it will attempt to grow through 
voxels of similar density first. 

Two-click (proximal and distal) vessel tracking with automatic centerline: 

There shall be a specific segmentation that attempts to identify the most probable 
centerline through a vessel given two points along that vessel. The user shall be 
able to hit a single button to enter the mode, it will prompt him to scroll to the 
appropriate slice proximal to the vessel section of interest and place a seed point 
using a mouse button click. The user will then be prompted to scroll to the 
appropriate slice distal to the vessel section of interest and place a second seed 
point using another click. The module will then attempt to find the centerline that 
most optimally connects the two points. 

Multiple-click vessel tracking: There shall be specific segmentation that utilizes 
more than two points to generate a vessel centerline that passes through all the 
points in succession. This is similar to the two-click tracking, but more user- 
intensive and more precise. 

Single-click vessel tracking: There shall be specific segmentation that utilizes 
exactly one point to generate a vessel centerline. It will find the longest vessel that 
passes through the initial point. It may not always find the vessel endpoints that 
the user has in mind, but it is useful for simple branch-less vessels. 
Vessel diameter and area measurements: The module shall make automatic 
calculations of vessel diameter and area over the length of any given centerline. 

o It shall display the vessel area or diameter (user selectable) in graphical 
form parallel to the main direction of the vessel. 

o It shall not include calcification in the area (or stenosis) measurement. 

■ The calcification shall be determined by a user-specified HU value. 
Percent stenosis calculation: The percent stenosis shall be calculated based on 
the area given a centerline, a normal region of the vessel, and a stenotic region of 
the vessel. 

o The user selects on the graph the location of normal and stenotic vessel 
and the percentage is then displayed above the stenotic region. 

o Taking a snapshot of the image allows capture of the calculation into the 
report. 

Centerline length calculation: There shall be the ability to measure the 
piecewise Euclidean length of a centerline either for the entire length or between 
two points on the centerline. 

Calcification Cleansing: There shall be the ability to exclude all calcification 
from the segmented vessel regions based on the user-selected HU threshold. 



• The module shall provide simplified selection of standard angiogram views 
including: 

o Left anterior descending views: 40° caudal 25° RAO, 40° caudal 40° LAO 

o Circumflex: 20° caudal 20° RAO, 0° ap 20° LAO 

o Right coronary: 30° caudal 30° RAO 

o Spider view for bifurcation/left main: 25° caudal 50° RAO 

• Generate reports specific to the patient with vascular measurements: 

o Evaluation reports. These reports shall provide data about the actual status 
of the arteries. 

• The views of the V3D Vascular shall include: 

o Selection View 
o Analysis View 
o Full-Screen Analysis View 

Selection View 

• Selection View shall include the following images: 

o 3D volume rendered image of the original intensities 
o 3D volume rendered image of the vessel-enhanced intensities 
o Orthogonal, multiplanar reformatted (MPR) images 
o 3D translucent heart view with coronary arteries depicted with removal of 
all ribs. 

o 3D view of pulmonary arteries with removal of all tissues that are not 
inside the lung region. 

• The selection view shall allow the following, in addition to the usual visualization 
and segmentation capabilities of the view: 

o Easy vessel tree selection from the 3D views 

o Creation of a new vessel from one, two, or more seedpoints 

o Creation of a new curved MPR centerline 

Analysis View 
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• The analysis view shall contain the following images: 

o A 3D overview image in the upper left 

o A curved MPR image on the right half. The curved MPR image can be 
any of the following types: 

■ Curved MPR 

■ Luminal MPR 

o A detail 2D/3D image in the lower left. It will always show a detailed 
view of the currently selected location along the vessel. It can be any of 
the following types: 

■ Oblique cross-sectional MPR 

■ Endoluminal 

■ Axial MPR 

■ Detail 3D (subcube centered at the current location) 
Full-Screen Analysis View 

• The full-screen analysis view shall contain a single full-screen image of one of the 
following types: 

o A 3D overview 

o A curved MPR image (any of the three types) 

o Oblique cross-sectional MPR 

o Axial MPR 

o Detail 3D subcube 



Planar Reformation of Vascular Central Axis Surface 

with Biconvex Slab 



Abstract 



Curved MPR is one of the popular vascular visualization methods. It re-samples and visualizes the vascu- 
lar central axis surface (VCAS) that is a curved surface passing through the vascular central axis (VCA) 
or vessel centerline. This results in a set of 2D images. In this paper we present a 3D method, VCAS pla- 
nar reformation (VPR) by a convex hull called biconvex slab. With a biconvex slab the entire vessel is 
enclosed and rendered in one image by volume rendering, such as MIP or X-Ray. The method is applied 
to CTA data sets. The resulting image is clear and free from obstruction of bones and other adjacent or- 
gans. In addition, since the biconvex slab minimizes the vessel volume, its rendering is in real time. 

Keywords: computerized tomographic angiography (CTA), vessel analysis, vessel visualization, curved 
multi planar reformation 

1 Introduction 

Five image reconstruction and visualization techniques have been widely used in Computerized Tomo- 
graphic Angiography (CTA): raw data axial, multi planar reconstruction (MPR), maximum intensity pro- 
jection (MIP), shaded-surface display, and volume rendering. Addis et al have compared these techniques 
[AHIL01]: volume rendering is an accurate method for evaluating all grades of stenosis. However, these 
methods are inadequate to visualize vascular structures. For instance, the entire vessel cannot be visual- 
ized in one image, including its lumen, wall, and surroundings. One way to visualize vascular structures 
is to resample and visualize the vascular central axis surface (VCAS) that is a curved surface passing 
through the vascular central axis (VCA) or vessel centerline. This process is referred to different terms as 
curved Multi Planar Reformation (curved MPR), Curved Planar Reformation (CPR), and Medial Axis 
Reformation (MAR). In the context of vascular visualization these terms are not precise enough to de- 
scribe the fact that VCA is located on this curved surface. Therefore, we use VCAS to identify this 
unique property: a curved cross-section that passes through the entire VCA. Hereby planar reformation is 
explicitly mentioned as a process to flatten the VCAS. By this technique - VCAS planar reformation 
(VPR), the entire vessel is flattened on a planar surface and the whole vascular centerline is displayed on 
a single image. 

Generally speaking VPR techniques allow the investigation of the vessel lumen in a longitudinal cross- 
section through the VCA. However, vascular abnormalities (stenosis and calcium) might not be scanned 
by this surface and therefore they do not appear in the generated image. One way to overcome this prob- 
lem is to rotate the VCAS along the longitudinal axis. This results in a set of 2D images. These 2D im- 
ages are used to diagnose calcification and stenosis as well as other vascular diseases. It works the same 
way as viewing 2D CTA slices to understand the 3D relationships and positions of objects. There is no 
3D information on the images. 

In order to visualize the 3D vessel in the VPR resulting image, a new vascular visualization method 
called VPR with biconvex slab is presented in this paper. Instead of resampling on a 2D thin surface, we 
investigate the VCAS in a thick biconvex slab; see illustration below. By introducing the biconvex slab, 
we can visualize the entire vascular volume in one image and minimize the obstructions from its adjacent 
organs, such as bones. 




VCA extraction is the basic procedure for vascular analysis. There exists a wide variety of VCA extrac- 
tion algorithms. Based on the input data we may categorize them into two groups, using segmentation 
data or using raw data. In the segmentation data group, there are maximum inscribed sphere [ShPB96], 
3D thinning algorithms based on the grass-fire definition [LeKC94], a minimum-cost path using 
Dijkstra's shortest path searching algorithm [BiKSOl] and inner Voronoi diagram [OgI192]. VCA extrac- 
tion using raw data is also called direct tracking, such as using Dijkstra's shortest path [FeWKOl], wave 
propagation tracking [QuKiOl], and intensity ridge method [AyBu02]. Kirbas and Quek have a review 
about vessel extraction techniques and algorithms [KiQu03]. In general VCA extraction algorithms find 
the vessel centerline and some other corresponding geometric information such as maximum and mini- 
mum diameters, contours, area, etc. at each point of the centerline. 




Traditional VPR: it is a 2D image. Biconvex Slab VPR: it is a 3D volume rendering. 

In this paper, we focus our attentions on how to render 3D information on VPR, which is how to visual- 
ize VCAS with biconvex slab. Section 2 describes related work on VPR. In Section 3 we introduce our 
biconvex slab in details. Experiments and comparisons are presented in Section 4. Section 5 concludes 
our method and proposes the future work. 

2 Related Work 

The research of VPR focuses on two points: 

• How to visualize entire vascular lumen and wall in one image; 

• How to visualize the entire vessel tree in one image. 

Ideally we would like to render the entire vascular lumen in one image. Kanitsar at al recently have pro- 
posed to use a helical scan line starting from center point to scan the vascular lumen instead of the 
straight scan line [KWFG03]. The resulting image of helical CPR rolls out the vascular lumen. It might 
visualize stenosis and calcification clearer than normal curved MPR. But it is really not easy to under- 
stand the 3D information from helical CPR, the right position and orientation of calcium and stenosis. 
The problem is caused by the 2D imaging of CPR. 
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Kanitsar at al [KaFW02] and He at al [HeDLOl] both have researched on how to put the entire vascular 
tree in one image. But their method still does not help the radiologist to find vascular abnormalities effi- 
ciently. There are many other clinical references of curved MPR in different vascular applications. 

The focus of this paper is on the first point, how to render 3D information in VPR to find out calcifica- 
tion and stenosis easily and precisely. 

3 Method 

Vascular central axis surface (VCAS) can be represented by a ruled surface in mathematics, i.e. a surface 
that can be swept out by a moving line in space, and has a parameterization form of, 

r(u,v) = a(u) + vl(u) 

where a(u) is a 3D curve called directrix (or called base curve) and l(u) is the director vector. The 
straight lines themselves are called rulings. 

For curved MPR, a(u) is the vascular central axis (VCA) and / (u) is a constant vector, i.e. the vector- 
of-interest (Voi). Voi is usually chosen as a vector orthogonal to the main orientation of the VCA. 
Thus, VCAS is rewritten as, 

VCAS(u,v) = VCA(u) + vVoi 

The Gaussian curvature of VCAS is everywhere zero, thus VCAS is developable. In another word, they 
can be flattened on a plane. The VCAS is filled by scanning and re-sampling each ruling in the volume 
data and thus creates a curved MPR. In order to view the entire vessel without overlapping, curved MPR 
may stretch VCAS along the main orientation of VCA (the longitude vector of the image) in different 
ways, such as stretched MPR, and straightened MPR. 

Traditional curved MPR is 2D imaging and lacks of 3D information of the entire vessel. To create a 3D 
VPR we need a slab, i.e. a thick VCAS. The straightforward idea is to create a thin slab by sweeping the 
VCAS along the view direction. But a thin slab has some disadvantages to render VPR. 

• A vessel is a thin object and is often located near other organs. A thin slab may include other 
adjacent organs. When the vessel has variant diameters, the thickness of the thin slab is very 
difficult to control. As a matter of fact, there are often too many obstructions that hide the 
views of the vascular lumen. 

• The vessel centerline is often very long. This results in a very long slab after stretching. Thus 
to render a very long slab becomes a time consuming task. 

A biconvex slab is designed to overcome these shortcomings and convey precise 3D spatial information. 
3.1 Biconvex Slab 

Basically we need a vascular central axis (VCA) and a vector-of-interest (Voi) to construct VCAS. At 
each point of the VCA, a straight line is defined by Voi; see the solid lines in the figure below. This is the 
scan line of VCAS to resample the volume. Apparently the results are highly dependent on the orienta- 
tion of Voi. For instance the scan line in the figure below misses both calcium. 
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In order to view the 3D vessel, we need to create a hull to enclose the entire vessel, which is named 'bi- 
convex slab' in this paper. On each cross-section that is through the center point and perpendicular to the 
centerline, a convex hull is created using the diameters; see the dash line in the image. The orientation of 
the convex hull is decided by Voi. If we connect all the convex hulls along the centerline, a biconvex slab 
is constructed. In order to create the convex hull, we need other geometric information, such as the 
maximum diameter at each center point. 

Since biconvex slab is a 3D volume, we may use volume rendering to view the entire vessel, such as 
MIP, X-Ray methods. Due to the fact that the resulting image of VPR is a flattened plane, only a parallel 
projection is considered for biconvex slab rendering. 

3.2 Image Space Based Convex Hull Construction 

Considering parallel projection, we can set up the image space of VPR by defining viewing vector as 
View = Up x Voi. Each scan line of the VPR image includes left point (L), center point (Q, right point 
(R), and maximum radius (r), where 

CR = Voi, CL = -Voi, | LR |= /e/ig* A(Scanline) 




Suppose the scan line LR is a thin ribbon, then we may rotate the strip 90 degree along Voi to view it on 
the plane of Voi and View. The length of vessel's projection on Voi is within 2r. Assuming that the orien- 
tation of the vessel's contour is unknown, the hull is defined as a bounding box that has a square shape 
size of 2r • 2r. Considering the margin 5 that transmits the slab thickness from 2r to thin ribbon, the scan 
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line LR is divided into three segments, of which LL H and R»R are the scanning range, and L H R H is the 
rendering range. For the scanning range, we calculate the value the same as normal curved MPR, i.e. re- 
sampling it in the volume, supposing its thickness is 1 voxel. 

For each pixel P located within the 3D rendering range, a ray is emitted through P along the View direc- 
tion. For the ray of which the distance to C (|CP|) is less than r, the rendering depth of the ray is within 
±r, i.e. (P-r • View, P+r- View). For the rays located in the margin, the depth is interpolated between r and 
1, supposing that the minimum thickness is 1 voxel. 

Since the purpose of VPR is to investigate the vessel lumen, the main volume rendering methods we use 
are MIP and X-Ray. There are two ways to flatten the slab, stretching the slab called curved VPR and 
stretching the centerline called luminal VPR. 

3.3 Minimum Convex Hull 

Square bounding box is a loose convex hull. There are two methods to minimize the convex hull, using 
volume data or using geometric data. 

To use volume data, i.e. vessel segmentation volume, the initial ray estimated by the square bounding box 
will traverse through the segmentation volume first to calculate the entry point (P entry ) and exit point 
(Pan)- Considering a margin 8, the final ray will accumulate the volume contribution within (P enny - 
8- View, Pexu+S- View). 

To use geometric data, such as contours or the orientations of maximum and minimum diameters, (if we 
have only orientations of maximum and minimum diameters, we may estimate a rough ellipse supposing 
that vessel contour is an ellipse), we project the contours along View direction to scan line (LR). A ID Z- 
buffer is used to the find both the maximum forward and backward depth (CQ-View) on scan line. 
Thus, each pixel on the scan line will have two depths, Df (plus - forward) and Db (minus - backward). 
Considering a margin S, the volume rendering region is (P - (Db+8) • View, P + (Df+8) • View). 



4 Experiments and Results 

We used clinical CTA data sets to demonstrate the efficiency of biconvex slab VPR compared to normal 
2D MPR and thin slab MPR. 

We tested our method first in a coronary artery CTA data set. The left anterior descending (LAD) artery 
is selected. LAD is rendered by a curved MPR with the normal 2D method, the thin slab MIP method 
and the biconvex slab MIP method. The maximum diameter of the LAD is about 5mm in this case. We 
put a 1mm margin on both sides. Thus, the thickness of the thin slab is about 7mm.The resulting images 
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are shown below. The 2D method cannot show all the calcification. Thin slab method displays the calci- 
fication, but it sacrifices the clear boundary of the vessel. Furthermore, some parts of the vessel are hid- 
den by its adjacent organs, especially at the top and bottom parts of the LAD. The biconvex slab renders 
all the calcification and keeps a clear boundary and overview. 
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From left to right (a) LAD curved MPR with normal 2D method, (b) 7mm thin slab MIP method, (c) 
VPR with biconvex slab MIP method 



The second test is an abdominal aorta CTA data set. Since the diameter of an aorta is much larger than an 
iliac artery, the maximum diameter is about 54mm in the aorta and about 10mm in the iliac in this case. 
In order to view the entire aorta and iliac arteries, we set the thickness of thin slab is about 58mm 
considering a 2mm margin on both sides. Compared to normal 2D method and thin slab method, 
biconvex method displays the entire vessel without any obstruction and blurriness from bone and other 
adjacent organs. 




Main Abdominal Aorta Curved MPR: From left to right (a) Normal 2D method, (b) flat slab MIP 
method, (c) VPN with biconvex slab MIP method 

In order to visualize the vessel lumen, we use the X-Ray model in volume rendering. In the MIP curved 
MPR, calcification obscures the lumen size, which is needed to diagnose stenosis, see images below. In 
order to show the lumen, electrical cleaning is used to remove the calcium. After calcium cleansing, we 
render the vessel volume with X-Ray model. In the resulting images, the real lumen size is clearly visual- 
ized in both curved VPR and luminal VPR. 
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From left to right (a) Biconvex slab MIP curved VPR, (b) Biconvex slab XRay curved VPR after 
calcium cleansing, (c) Biconvex slab XRay luminal VPR after calcium cleansing 



5 Conclusion and Future work 

In this paper we present a biconvex slab VPR method. Compared to normal curved MPR and thin slab 
VPR, we conclude from our experiments, 

• Biconvex slab supplies more 3D information than normal curved MPR. Calcium and stenosis are 
more easily identified relative to the 3D position and orientation. 

• Compared to thin slab VPR, biconvex slab visualizes vessels clearly and free from obstruction of 
adjacent organs. 

• Biconvex slab is rendered in real-time. VPR resulting image is usually bigger than a normal 3D 
rendering image, especially for a long vessel. If we use a thin slab, the image space volume ren- 
dering will be very time consuming. Biconvex slab efficiently focuses the volume-of-interest. 

We present the biconvex slab VPR for a single vessel. Our future work is to construct the biconvex slab 
for entire vessel tree and render the 3D vessel tree VPN. 
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Calcium Cleansing in CTA Vascular Visualization 



1 Introduction 

One of the advantages of CTA over MRA is its capability to view both calcification and stenosis. But in 
3D volume rendering, especially in MIP, calcium on the vessel walls always blocks the view of lumen 
since intensity of calcium is higher than lumen's intensity. This makes it difficult to visualize the vessel 
lumen, i.e. to view the stenosis in 3D, see the images below. We observe the calcium in both surface- 
based and MIP volume rendering images. But we have no idea the volume of calcium and the stenosis 
caused by the calcification. Vessel segmentation is to detect bright tubular objects on dark background. 
Thus, in most cases, calcium or parts of calcium are segmented as vessels. In MIP calcium always over- 
whelm lumen voxels in images. Therefore, we need calcium cleansing in CTA vascular visualization. 




Calcium in CTA Vascular Visualization 



2 Method 

In general, the intensity ranges of calcium and vessel lumen that is enhanced by contrast agents are dif- 
ferent. Calcium intensity is above 500. And lumen intensity is below 500. In different studies, they may 
vary due to the time control of contrast agents. Basically there is a threshold to separate the calcium and 
vessel lumen. 

It does not work well if we only remove calcium voxels with a threshholding. Due to partial volume ef- 
fects, the voxels near calcium have higher intensity than lumen. It still blocks our view to the lumen. We 
must filter the calcium region after applying threshholding. 

The solution is to remove the partial volume effects by a Gaussian filter. Supposing that each voxel is 
sampled within a region with a Gaussian mask, when a voxel is changed to another value, it affects its 
neighbors located within the sampling radius. Thus, when we fill the calcium regions with a background 
value, we should remove the partial volume effects by a difference caused by filling. The radius (3.0* 8) 
of Gaussian filter kernel is set to 5 = 1 .0 based on our experiences. 

This is the pipeline of calcium cleansing. 

1) Filling the calcium regions that is set by threshholding with a background value. 

2) Calculate the 3D Gaussian Volume Mask (G) 
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3) For each calcium voxel, 

• Calculate the difference (d) between background and original intensity 

• Calculate the difference of volume mask by d*G 

• For all voxels within the volume mask, if it is not a calcium voxel, remove the corre- 
sponding difference in difference of volume mask. 

3 Experiments and Results 

We use a carotid CTA data set to show the results of calcium cleansing. After calcium cleansing, a clear 
view of lumen is rendered in both images below. 




Vascular Visualization after Calcium Cleansing 



The 2D slices images show the results before and after the calcium cleaning. 




Calcium in CTA Vascular Visualization 



4 Conclusion 

In this paper we presetted a technique called calcium cleasning to remove the calcium within the vessel 
to clearly visualize the vessel lumen. The unique contribution is the removal of parital volume effects of 
caliums. This technique is very helpful for physicians to diagnose the stenosis in CTA. 
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V ssel S gmentation in CTA Data S ts 
with Vessel ness-based Front Propagation 



Abstract 



With the progress of multi-detector computerized tomography (MDCT) and increasing temporal and spa- 
tial resolution of data sets, clinical use of computerized tomographic angiography (CTA) is increasing. 
Vessel segmentation can be quite challenging, but is needed to isolate vascular structures. In this paper a 
fast and robust method for CTA vessel segmentation is presented using a two-step procedure, a pre- 
processing step and an interactive step. During the pre-processing step, a vesselness volume is computed 
by using a Hessian filter preceded by a CTA pre-filter. In order to enhance different size vessels, a MIP- 
volume pyramid is created and multi-level vesselness is computed and merged. Furthermore, the vessel- 
ness volume is integrated into a speed volume for front propagation based on vesselness, gradient and 
zero-crossing. At the interactive stage, a vessel central axis (VCA) is tracked and shown in real time by 
each click. Front propagation takes up to several seconds to finish segmentation after initialized by 
VCAs. The unique contribution of this method is the fast multi-level vesselness computation that com- 
putes the vesselness efficiently and the pre-CTA filtering that makes the vesselness applicable to CTA 
data sets. This method has been utilized in CTA vessel measurement due to its speed and accuracy. 

Keywords: vessel enhancement, vessel detection, volumetric vascular segmentation, level set. 



1 Introduction 

Vascular information is crucial in many medical applications such as in characterizing stenosis and aneu- 
rysm formation in cardiovascular diseases. The most commonly used medical imaging techniques that 
acquire vascular structures are ultrasound, computerized tomographic angiography (CTA), magnetic 
resonance angiography (MRA), and Digital Subtraction Angiography (DSA). In the past, CTA played a 
limited role in heart disease evaluation due to technical limitations, i.e. temporal resolution, necessary to 
glean essential anatomical, functional, and physiological cardiac data [BeOSOO]. With the introduction of 
multi-detector computerized tomography (MDCT), i.e. 16-slice MDCT, along with three-dimensional re- 
construction technology [MoHK03], it not only allows routine studies to be performed much faster than 
with single-detector CT scanners, but also makes available new applications, especially in the filed of 
CTA [HoFiOl]. CTA is gaining wider acceptance in clinic practice relative to the other modalities and 
has proved extremely valuable for many applications [FishOl], such as in cardiac imaging [Ache03, 
FlSk03]. 

Data acquisition with CTA requires proper timing of the contrast bolus but is otherwise straightforward 
and reproducible. Contrast enhanced (CE-) CTA that is essential a record of vessel lumen is analogous to 
CE-MRA but different from Gd-MRA [PYKH93]. Compared to MRA, CTA offers superior temporal and 
spatial resolution [WeinOl]. In addition, CTA can visualize various types of plaque including calcifica- 
tion and visualize the luminal blood through a contrast bolus. CTA images can also visualize the adjacent 
bony structures offering landmarks for surgical planning. MRA cannot be used for patients with pace- 
maker and certain metallic prostheses. Post-processing and data visualization play a significant positive 
role on CTA applications [FRRV01]. Due to bone and calcification, vessel segmentation in CTA, one of 
the most important post-processing, is much more complex than in MRA. It is the main challenge in ap- 
plying CTA in routine clinical use. 




Generally speaking, vessels segmentation in both CTA and MRA is to detect bright tubular objects in a 
dark background, e.g. the vessel centerline is explained as the ridge of intensity height surface [AyBu02]. 
But due to calcium, there are some high intensity spots (a type of hard plaque) within vessels in CTA 
data sets, particularly in elderly patients due to advanced atherosclerosis. Plaques are located within ves- 
sel walls and thus change the profile of local signal intensities. They can be mistaken as part of the vessel 
lumen (missing the real lumen) or as part of bones (missing the plaques). Besides, vessels may be very 
close or attached to bones at some body parts, such as the vertebral artery as it courses through the cervi- 
cal spine and the distal leg arteries in the region of the calf and ankle. In both situations (bone and cal- 
cium), vessel segmentation becomes very challenging. 

There is a large body of literature about vessel segmentation, including wave propagation [QuKiOl], 
multi-scale line filter [Lind94, SNSA98, FNHW99], deformable model [KaWT88, RBFM97], live-wire 
[FaUMOO], and level set method [LFGK01,MaSV95]. But very few papers discuss vessel segmentation 
in CTA cases. Most of the literature discusses MRA data sets. Kirbas and Quek have a review about ves- 
sel extraction techniques and algorithms [KiQu03], in which the primary application addressed is the 
segmentation of vascular structures in MRA and X-Ray Angiogram (XRA) data sets. Felkel has a review 
about vessel segmentation of peripheral vessels in CTA [FelkOO]. CTA has its own context different from 
other imaging modalities. Recently due to the widening acceptance of CTA in vessel analysis, vessel 
segmentation in CTA draws more and more attentions from researchers. Our research in this paper fo- 
cuses on Hessian filter and front propagation applicable for CTA. 

Multi-scale Hessian filter has been proposed to enhance tubular structures in three-dimensional medical 
images [SNSA98, FNVV98]. This method involves convolving the image with Gaussian filters at multi- 
ple scales and analyzing the eigenvalues of the Hessian matrix at each voxel in the image to determine 
the local shape of the structures in the image, such as tubular structure, planar structure, speckle noise, 
and other structures. Some specific functions are designed with the three eigenvalues and eigenvectors of 
Hessian matrix to enhance the tubular structures, so called "vesselness". This vesselness volume is di- 
rectly used to visualize and detect vessels. To detect the intensity ridge with Hessian matrix [AyBu02] is 
another way to use multi-scale Hessian filter to track vessel central axis and estimate its width. Both 
methods use the same mathematic theory, i.e. Hessian matrix. 

The front propagation method [MaSV95] uses the level set method developed by Osher and Sethian 
[OsSe88] and adapts it to three-dimensional vascular segmentation. The propagating interface is con- 
trolled under a curvature dependent speed function. The level set method can easily handle the changes in 
topologies and is independent of the parameterization of the evolving surface model. 

We choose CTA data sets as the context of our vessel segmentation. The paper proposes a fast vesselness 
computation method adaptable for CTA datasets. Based on vesselness, we present a fast front propaga- 
tion method to segment vessel with minimal user interaction. Our method shows promise in CTA vessel 
measurement due to its speed and accuracy. The subsequent sections focus on the following aspects of 
our method: 

• Section 2 discusses a fast multi-level vesselness computation method. It discusses the principal 
difficulties in vesselness computation and introduces our multi-level vesselness computation 
method by using a MlP-volume pyramid. In addition, it deals with how to get the maximum re- 
sponse by filtering CTA data sets, called CTA-prefiltering, before applying Hessian filtering. 

• Section 3 presents a fast front propagation method based on vesselness by using a time-tag nar- 
row band. It discusses how to construct the speed volume based on vesselness and how to mini- 
mize the user interaction by using vessel central axis (VCA) tracked in speed volume. 

• Finally in Section 4, we present our experiments and results. Five CTA studies in different parts 
of the body are selected and tested. Experimental results including run time are presented. It 




takes several seconds to segment a vessel. The results are accurate especially in situations where 
traditional methods have difficulty in separating bones and blood vessels. 



2 Computation of Multilevel Vesselness 



Generally speaking, a vessel is considered as a linear or tubular structure. Using a multi-scale line filter 
to detect or enhance tubular structures is one of the most popular methods, especially in 3D. Vesselness, 
named by Frangi [FNVV98] and Sato [SNSA98] in their own papers, is a grade to assess tubular struc- 
tures. The higher vesselness a voxel has, the more probability it belongs to a tubular structure. Although 
not all of the tubular structures in vascular imaging are vessels, we can discriminate vessels from their 
soundings with vesselness. 

Supposing a vessel model is a cylinder in which the Z-axis is the ves- 
sel central axis and x-y plane is the vessel cross-section with Gaus- 
sian distribution, I 0 is that image, see Figure 1 . 
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Eq. 1 



The Hessian matrix H is given by 
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The three eigenvalues and eigenvectors of the Hessian matrix //are 
4 =0,7, =(0,0,1); Eq.2 




Figure 1 Vessel model with a 
Gaussian luminance profile 



K = ~ $ , v 2 = (x, y,0) ■A 3 =--±,v 3 = (-y, x,0) ; 

in which V,, V 2 and V 3 are vectors perpendicular to each other. 

Based on this model, vesselness is scaled by the linearity of ^ with X, and A,, i.e. /l^Oand 
Aj < %2 « 0 , and the symmetry of and , i.e. ^//^ « 1 . 

Different functions may be designated for S^e and 5 symm . Frangi used exponent function [FNVV98] and 
Sato used polynomial function [SNSA98]. 

The main issue of vesselness computation is that of computational cost. The multi-scale convolutions of 
2 nd derivatives can be very computational costly. For example, Sato took 10 minutes to calculate a three- 
scales vesselness in a MRA volume size of 256x256x102 with eight CPUs 168MHz Sun UltraSparc 
processor, regardless of the thresholds used to remove the background. Frangi spent similar time in their 
paper with only a single scale filter. How to reduce the computational cost without loss of vesselness 
quality is one of the main challenges in vesselness computation. 

On the other hand, vesselness is only mentioned as a vessel enhancement method used in MRA data sets 
[KiQu03, FelkOO]. Both authors, Frangi and Sato, experimented with their methods using MRA data sets. 
Sato used a head MRA case and Frangi used a Carotid MRA case. In a later paper from Sato [SWBN00], 
a bronchial case with threshold -180HU and a liver CTA case with threshold range of 0HU ~ 300HU 
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were discussed. Both cases had no adjacent bony structures and ignored calcium. Since some bones have 
tubular shapes, such as ribs compared to aorta, how to reduce the "false-positive" vesselness is another 
challenge to computing vesselness in CTA. 

2.1 CTA Pre-filtering 

In CTA scanning, an iodinated contrast bolus is injected into a large vein and rapidly circulates through 
the heart and then reaches the arterial vessels. The histograms below show the effect of an IV contrast 
bolus in an abdominal CTA case. We show the histogram range from -100HU to 900HU in Figure 2. The 
dash line is the histogram before injecting contrast agent. The solid line is the result of contrast enhanced 
scanning. From the overview of both histograms, we can observe the shift of voxels in a low intensity 
range [-100, 0] to a high intensity range [50,500], see also the detail view in Figure 2. That is to say, this 
range of Houndsfield units is markedly changed (enhanced). 




Figure 2 Comparison of histograms from an aorta CTA scanning with- and without- contrast 
agent. Dash line is the histogram scanned before injecting contrast agent. Solid line the histo- 
gram scanned after injection. 

An iodine contrast agent does help to distinguish vessels, but a vessel's intensity range still overlaps with 
other structures, especially lower-density bone and marrow. For example when the vertebral artery goes 
through the cervical spines or anterior tibial artery is in close proximity to the tibia, they have a similar 
intensity range as bone surfaces (low-density bone). Due to timing control and circulation (caused by 
plaques in vessels) of the contrast agent, the intensity of enhanced vessels is quite variable, from 100HU 
to 600HU. Because of this, a simple threshholding will not work well or consistently in CTA cases. 

The other obstacle of CTA vesselness computation is that in heavily diseased arteries, calcium and other 
hard plaque adheres to the vessel wall and changes the local intensity profile, which makes it very diffi- 
cult to compute the right eigenvalues and eigenvectors of the local Hessian matrix. Because high-density 
bone and calcification share the same intensity range, vesselness becomes discontinued or broken as a re- 
sult of threshholding. Ideally, we would like to keep the calcium within vessels as parts of vessels seg- 
mented. But a Hessian matrix will have difficulties in getting the right response when encountering hard 
plaques. 

As we mentioned at the beginning of this section, the ideal vessel model is a Gaussian luminance profile, 
see Eq. 1, i.e. the highest intensity is at the middle of the tube and the intensity declined based on Gaus- 
sian function at the boundary. But this does not occur in CTA clinic practice, and we cannot directly ap- 
ply the Hessian filter. Therefore, CTA data sets need a pre-filtering step to enhance the potential tubular 
structures before calculating vesselness. This CTA pre-filter should satisfy the following requirements: 



• Keep the Gaussian shape vessel luminal profile. 

• Adjust the volume intensity so that the maximum intensity within the vessels lumen becomes the 
maximum intensity of the volume. 

• Normalize the intensity in order to compare the vesselness from different locations. 

The CTA histogram can be categorized into three ranges: Ex-vessel Low (ExL) range including air, fat 
and soft tissue, In-vessel (In) range including contrast enhanced vessel, low intensity bone and marrow 
etc., and Ex-vessel High (ExH) range including bone and calcium (see Figure 3 below). We observe that 
at both ends of each range there is a Significant Point (SP). At this point, the slope changes significantly 
because the Region In Range (RIR) reduces considerably. SP is regarded as the statistical threshold to 
sort the different materials in histogram. The SP is calculated by the intersection of two asymptotes (ap- 
proached by tangent lines). In each range we can calculate a pair of SP points A(SP 0 , SPJ. From SP 0 the 
corresponding RIR starts growing massively. After SP t the growing of RIR vanishes. 

Since the histogram is separated into three ranges, three SP point pairs are calculated, i.e. A EliL (SP 0 , SPJ, 
Am (SP Q , SPx), A Ex h(SP 0 , SPO-The pre-filter is set up as a roof-shape curve, specifically 

• Set SPi in A ln at the Peak Point 

• Set SP 0 in A ln at Left Verge Point 

• Set SP 0 in A Efii at Right Verge Point 

A quadratic curve, called the Normalizing Roof Curve (NRC), is approximated by these three feature 
points. A look-up-table based on NRC is calculated and used to map the original volume to keep the In- 
range voxels and dehance the ExL-range and ExH-range voxels. 



Histogram 




Figure 3 Three CTA intensity ranges and the Significant Point in each range 



2.2 Multi-level Filtering 

The linear filter is basically a multi-scale filter, i.e. it filters with different radii (cr) convolved at each 
voxel in the volume and the maximum response is chosen as the vesselness at the current voxel. But gen- 
erally speaking, a vessel's radius can differ by as few as several voxels (coronary artery) and up to a hun- 



dred voxels (abdominal aorta). Therefore, a multi-scale filter must calculate all sizes in order to find the 
maximum response. On the other hand, large-scale filtering is very time consuming. 

In a multi-level filter, a volume pyramid is created based on the original volume, called the MIP (Multum 
In Parvo) Volume. A MIP structure is widely used in computer graphics for 2D texture mapping, (MIP 
map) [Will83]. In a MIP Volume, Level 0 is the original volume. The next level volume stores voxels re- 
cursively with half of the pre-level volume size, i.e. filtering or averaging over every 2x2x2 voxel. This 
process can reduce the volume size to one voxel or a certain level. With this volume pyramid we can re- 
construct any volume of which its size locates between two integral levels. Every voxel can be interpo- 
lated with 16 nearest voxels to its neighbor integer level volumes (8 voxels on each level). 
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Instead of using a multi-scale filter, a filter with fixed scale is used in multi-level filtering. The size of the 
filter is set up based on the smallest vessel expected in the original volume, for example 5 -voxel -diameter 
vessel. Then the same filter is applied to compute the vesselness in other levels, such as 1.5, 2.0, ... etc. 
The resulting volume is scaled down to the 0-level size. Finally, results from all levels are merged into 
the 0-level volume by maximum operator. 

V = max(V,), / = 0,0.5, 1.0, 1.5, Eq.3 



Calculating the vesselness using a multi-level filter avoids the time-consuming convolution of large-scale 
filters. The disadvantage is that it may introduce blurring in the resulting vesselness volume. As we will 
see in the next section, blurriness within vessels won't affect propagation or segmentation. If we can 
keep a sharp vesselness near boundary, segmentation will get the same results as a multi-scale filter. In 
order to reduce the blurring near boundary we use a gradient weight function before high-level vesselness 
is scaled down to 0-level. 

v,° =(1.0- gradient) * v, Eq. 4 



3 Front Propagation based on Vesselness 

Front propagation [MaSV95] is an efficient fast marching level set method to track the monotonic ad- 
vance of interfaces with a speed that does not change its sign, either expanding or shrinking. In this sec- 
tion we introduce a front propagation in which the speed function is defined on vesselness. 

Briefly we can formalize the front propagation as the evolving function y/, namely 

y/, +F|V^| = 0 Eq.5 



where F is the speed function. 



Let F is defined as the implicit function such that F, ={ ^, + F\ V = 0 } , it is easily shown 



^ = -F\ V u/\ , where — = -FN, N = — 
dt 1 n ' dt V 



Eq. 6 



Thus, the evolution of zero level-set of y/, equals the evolution of F„ a time-dependent implicit surface 
representing the evolving segmentation interface. 

Considering the special case of monotonically advancing surface, i.e. front propagation, speed function F 
is always positive and normal vector N is always pointing outwards. Instead of numerical approximating 
the derivatives in Eq.5, we explicit construct and update the front interface F t with Eq.6. 

3.1 Timer-tag Narrow Band 

The point to explicit construction of F t is to create an active "working zone", called Narrow Band 
[MaSV95], a local region around the front. The narrow band is constructed by active seeds. Each seed is 
defined by a vector of (P, t), of which P is the target position and t is the Timer. It explains that the cur- 
rent seed will take f-time to arrive target position P. When t is positive, this seed is active. When t is less 
or equal than zero, this seed is deactivated and merged into r t . 

The narrow band is maintained by a heap data structure sorted by timer t, called front-seed heap. The 
seed with smallest t is at the top of the heap. At each time step, we check the heap and the front interface 
is marching outwards as below 

1 . For all seeds, decrease t 

2. For all the seeds of which t < 0: 

• remove it from heap 

• merge into r t 

• insert its neighbors into the heap 

3. Check the left seeds in heap of which t > 0: 

• If P is in F t , remove the seed from the heap 

For each neighbor point NB of front point P, such as 6-neighbors, 1 8-neighbors, or 26-neighbors, we 
check only the outside NB, i.e. outside the f, and initialize its timer with Eq. 7. 



We store the NB position and t in the seed and insert it into the front-seed heap. 
3.2 Speed Function 

Speed function defines the motion of front. In order to segment vessels, we need to design a proper speed 
function, which increases closing to the central part of vessels and vanishes at the boundary of vessels. 
Most speed functions are defined directly on original images. Here we introduce a speed function defined 
based on vesselness. 

Two motions affect the speed function F in Eq. 7, motion by curvature term K v and motion by the image 
term F/. 



t = — !—± — , where L = NB - P , F is the speed at point P 
(L-N)F 



Eq. 7 




F = {\.0-eK v )F,, 



Eq.8 



where K v is mean curvature and AT = V • ( — ) , e is the front smoothness control factor. 

I Vy\ 

Mostly F, term is directly calculated from original image by gradients, called static speed. Considering 
the vesselness, gradient and zero-crossing, a speed volume is initialized by using these two volumes. 

F — (I 

' 1 0 + C | grcid | vesie ' nesj com ) Eq. 9 



where 7 cons is a constant term to keep the front moving. Near the zero-crossing points, this term is van- 
ished. 

Since gradient and zero-crossing points are calculated while computing vesselness, instead of calculating 
F, on fly, a speed volume is integrated and saved directly after vesselness is calculated. By using the gra- 
dient and zero-crossing, we convert the vesselness that might be blurred by the multi-level filter to a 
sharpening speed at object boundaries. 

3.3 Front Initialization 

The front is initialized by the vessel central axis (VCA) tracked in speed volume. In order to minimize 
user interaction, VCA is tracked based on the primary direction V 0 in Hessian matrix (Eq. 2) from a use 
click. The VCAs along both primary directions (±V 0 ) are tracked. Unlike most ridge-tracking algorithms 
[AyBu02], we use single-scale Hessian filter to estimate the primary direction, the same scale as we used 
to calculate the vesselness. Advantages include: 

• The vesselness or speed volume is already "centralized" single-scale is usually enough to track a 
VCA. 

• VCA is used as initial front to segment vessel. Basically if the initial front is located within a 
vessel, the results segmented by front propagation keep the same. That is to say, initial position is 
not critical to segmentation results. 

In addition, the VCA can be trimmed off when two central axes are close enough in space, i.e. merged 
into one central axis. The VCA always follows the directions of user clicks. At bifurcations, the user can 
control the vessel segmentation by selecting the interested branches. Since the VCA is tracked and dis- 
played in real time, it gives the user an overview to interrogate the vessel that is going to be segmented. 

4 Experiments and Discussion 

The method is integrated into the vascular measurement system. The experimental data sets in this sec- 
tion are selected from our validation data sets collected from different clinical institutions and acquired 
by different CT scanners. CTA scanning from different parts of body are chosen, including the Head, 
Neck, Heart, Abdomen, as well as Lower Extremities (runoff). We demonstrate our method by segment- 
ing small radius vessels (Coronary Artery), large radius vessels (Abdominal Aorta) and long vessels (pe- 
ripheral arteries). In addition, we show the segmentation of the vertebral artery and the anterior tibial ar- 
tery, vessels that are in close proximity to bone. 

4.1 Experiments Pipeline 

Before vessel segmentation, all data sets are put through a five-step pre-processing. 



1) The original volume is interpolated into an isotropic volume. 

2) The isotropic volume is filtered by the CTA pre-filtering. 

3) The MDP-Volume is created and Hessian filter is applied to different levels. 

4) The responses from different levels are merged into the resulting 0-level volume. 

5) The results from the Hessian filter are integrated into the speed Volume 

The computation of second derivatives of Gaussian in step 3) was implemented by three separate one- 
dimensional convolutions. The radius (3 <5 f ) of Gaussian filter kernel is set with df=l.O. The pre- 
processing is launched automatically when the system server receives CTA data set. Based on the body 
part of the study identified in the DICOM header, the system selects the vesselness level and the voxel 
units of isotropic volume from a presets library. 

Table 1 lists the computational cost of pre-processing, including the size of the original and isotropic data 
sets and the total time of pre-processing. The raw data slices are 512x512. The computational computer 
is a 2.8GHz 2GB memory PC. Most of the pre-processing times are < 5min. Only the peripheral case 
takes about 14mins due to the large number of slices (about 1400 slices). 



Table 1 Pre-processing computational cost of experiment data sets 



Data Set 


Original Volume 
Data Sets 


Isotropic Volume 
Data sets 


Pre-processing Time 


Circle-Of- 
Wills 


512 2 x68 
(0.43 2 x0.75mm) 


512 2 xll7 
(0.43 2 x0.43mm) 


2m03s 


Carotid 


512 2 x336 
(0.488 2 x 0.50mm) 


500 2 x336 
(0.50 2 x0.50mm) 


5m59s 


Coronary Ar- 
tery 


512 2 x256 
(0.488 2 x0.40mm) 


512 2 x209 
(0.488 2 x0.488mm) 


4m58s 


Abdominal 
Main Aorta 


512 2 xl40 
(0.7148 2 x 3.0mm) 


333 2 x380 
(l.l 2 xl.lmm) 


3m34s 


Peripheral Ar- 
tery 


512 2 xll34 
(0.65 l 2 x 1.0mm) 


512 2 xl403 
(0.651 2 x0.651mm) 


13m39s 



In our vascular system, vessel segmentation is a part of the interactive vessel measurement process. It 
starts by a user click to select a vessel. After one click, one VCA is tracked and shown in real time. If 
there is more than one VCA, the VCA voxel searches other VCAs in its 26-neighbors to merge or trim 
off VCAs. If user thinks the shown VCAs adequately describes the vessels to be segmented (or to meas- 
ure), then he/she starts the segmentation just by pressing 'Enter'. The segmentation as well as the meas- 
urement finishes in several seconds. We investigate all selected cases one by one in the rest of this sec- 
tion. 

4.2 Circle-Of-Wills (COW) from Head CTA 

In this experiment, three levels vesselness are calculated, i.e. /=0, 1, 2. The results from all levels are 
merged into the 0-level volume; see the left image in Figure 4. Vesselness volume enhances most the 
vessel structures in the CT volume. 

By clicking on one of the vessels, a VCA is tracked immediately by the primary direction in Hessian fil- 
ter; see the middle image in Figure 4. It takes about 0.3sec to find the 275 voxel long VCA. The small 



23 



sphere in the middle of VCA is the point clicked by user. 'F' and 'B' are two ends of the VCA. The user 
may edit two ends by deleting voxels on VCA in case it tracks longer than expected. Front propagation 
takes about 3.2sec to segment the vessel in the right image in Figure 4. 



4.3 Carotid and Vertebral Arteries from Neck CTA 

In this experiment, the vesselness level considered is /=0, 1, 2; result sees Figure 5. After a single click it 
took about 5.8sec to segment the left carotid artery, marked as red in the resulting image. Notice that cal- 
cium in the carotid wall is included as part of the segmentation result. By a few clicks (in this case, 4) the 
VCAs of a vertebral artery are tracked. We can see the merges of VCAs in the image. It took 3.2sec to 
finish the segmentation. Since the vertebral arteries traverse the vertebral canal within the cervical spine, 
they have very similar intensity to the adjacent low-density bone. The slices show the details of the seg- 
mentation result. It is very difficult to segment them just based on intensity. Vesselness assists us in sepa- 
rating them based on their local shape, i.e. tubular structure. In addition, we show both left and right ca- 
rotid and vertebral arteries segmented. 



4.4 Coronary Arteries from Heart CTA 

In this experiment, the vesselness level is computed until 1, i.e. 1-0, 1. By a single click we find the ini- 
tial VCA of LAD artery in 0.43sec, which is 388 voxel long. The segmentation took 2.6sec; see the re- 
sults in Figure 6. It segments LAD artery from the aorta to its periphery. It demonstrates the ability of 
vesselness to find small radius vessels like the LAD artery. Due to the circulation, the peripheral portions 
of the LAD artery are poorly enhanced. On the slab MIP in the Sagittal view, the end of LAD artery is 
marked with translucent red. The intensity is very close to soft tissue nearby. According to our experi- 
ences, the intensity there is usually around 50HU. Again, vesselness helps us to separate it from sur- 
rounding tissues. With several clicks, the left arterial tree is easily segmented. 




Vesselness Volume 



Vessel Selection 



Vessel Segmented 



Figure 4 Circle-Of-Wills (COW) Experiment 





Figure 5 Carotid and Vertebral Artery Experiment 



Upper (from left to right): Vesselness volume, Carotid segmented, VCAs of vertebral artery, Vetebral 
artery segmented. 

Lower (from left to right): four cross sections of the vertebral artery segmentation result; carotid and 
vertebral arteries segmented 




Figure 6 Coronary Artery Experiment 



From left to right: Vesselness Volume and the VCA of LAD; LAD segmented, Segmented LAD in Sag- 
ittal MIP; Left Coronary Arterial tree segmented; 
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4.5 Abdominal Aorta and Iliac Art ries from Abdominal CTA 

In this experiment, the vesselness level is computed until level 3, i.e. /=0, 1, 2, 3. With this experiment 
we demonstrate the vesselness on big radius vessels. By two clicks, one on the aorta and the other on the 
left iliac, we created two VCAs. The second one from left iliac merges into the first one. The first VCA 
is tracked in 0.46sec, which was 362 voxels in length. The second VCA took about 0.3sec, which was 
257 voxels in length. The segmentation took about 8.1 sec. The results are shown in Figure 7. In the mid- 
dle of aorta, there are some white spots, which is the calcium. For cross sections with calcified walls, 
there are two contours. The red contour is the segmented vessel region and the green contour is the vessel 
lumen region 
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Vesselness Volume Aorta Segmented One Slice of Result 

Figure 7 Abdominal Aorta Experiment 



4.6 Peripheral Arteries from Lower Extremity CTA 

In this experiment, due to the large number of slices (1403 slices) and the limited memory available in a 
standard Windows (2.0GB) operating computer, we split the volume into two slabs and calculated the 
two slab's vesselness separately. At end the two vesselness slabs are merged into one vesselness volume. 
This experiment is designed to show the vesselness computation for extra long vessels and large vol- 
umes, especially for the full body CT scanning. 

By a single click in the vesselness volume, the VCA of right peripheral artery is easily tracked. It took 
1.4sec and the length is 1350 voxels. The segmentation took about 8.6sec. Resulting images see Figure 8. 
In addition, we show the whole peripheral artery tree. 
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Figure 8 Peripheral Lower Extremities Experiment 

Left: Vesselness volume and the one click VCA of right peripheral artery 
Middle: The segmented right peripheral artery 
Right: All peripheral arteries segmented 
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Another experiment in this data set is to show the segmentation of the anterior tibial artery, see Figure 9. 
Because the anterior tibial artery runs along the tibia and sometimes is contiguous to the bone, both in- 
tensity (vessel and bone) are very close after vessel enhancement. The detail view displays the situation. 
With normal region growing it is very difficult to segment it from bone. But vesselness helps us to sepa- 
rate vessels from the bone based on the local tubular shape. The result both on the 2D slice and 3D ren- 
dering image show a successful segmentation. 




Figure 9 Tibial Artery Experiment 



Upper left: the segmented Tibial artery 
Upper right: the detail before segmentation 
Lower left: a 2D slice of segmented tibial artery 
Lower right: the detail after segmentation 



5 Conclusion 

In this paper we present a vessel segmentation method applicable for CTA data sets. From our experi- 
ments, we conclude: 

• Creating vesselness with a Hessian matrix is a powerful tool to enhance the tubular structures in 
the volume. In CTA, since the intensity of enhanced vessels overlaps the intensity of low-density 
bone, vesselness is an efficient tool to separate vessels from nearby bony structures. 

• Calcium is often deposited in arterial walls. Segmented vessels should include both lumen and 
calcium. Otherwise, segmented vessels will be fragmented. 
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• Bone or calcium changes the local profile of the signal intensity, which reduces the responses of 
vesselness from a Hessian matrix. In order to get a better response, we use histogram analysis to 
design a Normalizing Roof Curve to pre-filter CTA data preceding the Hessian filter. 

• Instead of using a multi-scale filter, we use a multi-level filter and a MlP-Volume pyramid struc- 
ture to compute vesselness. Large size vessels are enhanced in high levels and their vesselness is 
merged, resulting in a 0-level volume. Hence we avoid the time-consuming convolution of a 
large radius filter. The computation cost in our experiments is reduced to about 5min in a stan- 
dard windows computer. 

• A unique speed function defined with vesselness is designed for the front propagation to segment 
vessels interactively. In order to minimize user interaction, VCA is tracked and shown in real 
time. By merging and trimming, users can directly view the vessels that are going to be seg- 
mented. The front interface is initialized by VCAs. Vessel segmentation is finished in a matter of 
seconds. 
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Optimization of Vessel Centerlines 



Abstract 

Vessel quantification needs an accurate and reproducible vessel centerline. Currently, most vessel center- 
line extraction algorithms are not stable enough for clinical applications. This paper presents a progres- 
sive optimization algorithm to refine a centerline after it is extracted. A new centerline definition for re- 
finement is proposed in terms of the minimum cross-sectional area. A centerline is divided into a number 
of segments. Each segment is corresponded to a local general cylinder. A reference frame (cross-section) 
is set up at the center point of each local cylinder. The position and the orientation of the cross-section 
are optimized within each local cylinder by finding the minimum cross-sectional area. After all center 
points are optimized, a refined centerline is approximated and re-sampled to the refined set of center 
points. This refinement iteration, local optimization plus global approximation, converges to the optimal 
centerline, yielding a smooth and accurate curve central axis. 

Keywords: medical visualization, vessel quantification, vessel skeletonization, medial axis transforma- 
tion, centerline 
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1 Introduction 

Analysis of vascular structures acquired by computerized tomographic angiography (CTA) or magnetic 
resonance angiography (MRA) is commonly performed for clinical diagnosis of vascular disease, e.g. as- 
sessing and monitoring stenosis secondary to atherosclerosis, for surgery planning, etc. Vessels can be 
evaluated using CT and MRI quantitatively - stenosis can be calculated by ratios of minimum to normal- 
ized diameter or cross-sectional area. Blood vessels can also be evaluated qualitatively using volume and 
surface rendering post-processing. Based on the tubular shape of vessels, a geometric model for vascular 
quantification utilizes a centerline and a series of cross-sections perpendicular to the centerline. Cross- 
sectional diameters and areas can then be calculated. An automatic reproducible vascular quantification re- 
lies on an automatic, reproducible and accurate centerline. 

The process to extract vessel centerline and its associated cross-sections is called vessel skeletonization. 
The term skeletonization and its generating Medial Axis Transformation (MAT) in 2D was first introduced 
by Blum [5]. Skeletonization simplifies the shape to the closest set of centers of maximal inscribed disks, 
which can fit within the object. The central locus of the centers is made the centerline. A 2D centerline 
with the radius of its associated disk allows an error free reconstruction of its tubular structure. Consider- 
ing the maximal inscribed sphere and the same MAT definition in 3D, a medial surface rather than a me- 
dial line can be found in objects [6], see Figure 1. The definition of centerline in 3D can then be ambigu- 
ous. Different definitions of skeletonization in 3D have been developed to find one line on this medial sur- 
face to be the centerline in 3D, such as maximum inscribed sphere [25], flow-based model using grass-fire 
transformation [15], and inner Voronoi diagram [20]. An interesting overview on different definitions of 
skeletonization using paradigms from geometry and mechanics can be found in [12]. 




There exists a wide variety of 3D skeletonization algorithms based on different definitions and extraction 
approaches. In the context of vessel skeletonization, many centerline extraction methods have been devel- 
oped. There are three basic approaches to centerline extraction based on input data: binary data, distance 
map and raw data. A good skeletonization [19] preserves the topology of the original shape, and approxi- 
mates the central axis. It is thin, smooth and continuous, and allows full object recovery. 

3D thinning algorithms based on the grass-fire definition of skeletonization are the main approaches to 
centerline extraction using segmented binary data. The boundary of a binary object is diminished itera- 
tively by deleting points that fulfill certain geometric and topologic constraints in each iteration step. Digi- 
tal topology and mathematical morphology are the basic theory applied in this thinning process. Topologi- 
cal thinning has been applied to vessel extractions by several authors [14, 16, 21, 9, 17]. Due to the sym- 
metric thinning, thinning processing usually does not reach a non-symmetric center point. Usually it con- 
structs voxel-wide graphs, which is not suitable for geometrical quantification. 

Centerline extraction based on a distance map is a direct application of Blum's skeletonization definition 
on binary data. The distance transformation of binary data creates a distance map, in which each voxel 
within objects has the minimum Euclidean distance to boundary [8]. The centerline is defined as the set of 
centers of maximal inscribed spheres constructed with the distance as its radius. A good survey on Euclid- 
ean distance transformation and its influence on vessel skeletonization can be found in [22, 19]. The proc- 
ess of searching or connecting a centerline computes a minimum-cost path in the distance map using 
Dijkstra's shortest path searching algorithm [3, 18, 6, 29]. In a 3D volume, voxels are taken as the nodes in 
a graph and neighboring voxels are connected by links. Each link is associated to the minimum distance to 
its boundary. The more centrally located a voxel; the lower its weight (defined as cost in a graph). The 
shortest (minimal cost) path is found and centered to get the vessel centerline. 

Direct tracking in raw data is another widely used technique [1, 2, 10, 23, 26, 27, 28]. Fellkel et al. [10] 
use Dijkstra's shortest path to track a path between two end points. The cost of links is associated with the 
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voxel value, e.g. absolute difference of the node value. Wave propagation tracks the normal path of the 
waveform from a seed point to detect the vessel centerline [23, 28]. Intensity ridge method [1, 2] converts 
an N-dimensional image to a surface in (N+l)-Dimensional space by mapping the intensity to height. The 
vessel centerline can be extracted via traversal of the one-dimensional height ridges on its surface. Frangi 
et al. [11] create a vessel geometric model using the snake model and refines it by optimizing the energy 
function. The central axis of the snake is the vessel centerline after refinement. 

Critical to any vessel centerline extraction technique is how well it handles noisy data, branches, and com- 
plex blood vessel anatomy. Generally speaking, centerline algorithms detect bright objects on dark back- 
ground, e.g. the centerline is explained as the ridge of intensity height surface. But due to calcification, 
there are some high intensity spots (known as plaques) within vessels in CTA data sets, particularly in eld- 
erly patients due to advanced atherosclerosis. Plaques are located within vessel walls and thus change the 
profile of local signal intensities. They can be mistaken as part of the vessel lumen (missing the real lu- 
men) or as part of bones (missing the plaques). A centerline should be centered based on the vessel walls 
and should also not break or twist due to obstructions caused by plaques and/or high-grade stenoses. Most 
of the current centerline algorithms have difficulties overcoming plaques in CTA studies. 

Another reason that the normal, discrete one-voxel-wide (some half-voxel-wide) centerline is not satisfac- 
tory in clinical applications is the non-reproducibility of vessel quantification. Quantification relies on an 
accurate and reproducible centerline. In fact, when one vessel is measured by different users or measured 
at different times or measured by different algorithms, the centerline may vary. This non-reproducibility or 
inaccuracy of quantification weakens its clinical application. Hence, in order to attain reproducible quanti- 
fication, centerlines need to be optimized to approximate the central axes, i.e., a good skeletonization [19]. 
Most current algorithms use smoothing after centerline extraction in order to remove the jagged changes in 
the centerline. But smoothing does not maintain centralization of the vessel skeleton (see one example in 
Section 5) in extracting the true centerlines in CTA studies. 



Wink and Niessen [26] discuss an algorithm to find a centerline by tracking the center of vessel contour 
on each cross-section, of which the orientation is calculated by its preceding points. In fact, the preceding 
and subsequent center points both affect the cross-sectional orientation. For a uniform cylinder, orienta- 
tion does not greatly effect the position of a center point, as shown in Figure 2. Point A is the center of 
cross-section AO that is perpendicular to the centerline. Although neither of cross-sections Al and A2 is 
perpendicular to the centerline, point A is still the center of both Al and A2 due to the symmetry of uni- 
form cylinder. For non-uniform cylinder, orientation is critical to the position of a center point. In Figure 
2 point B is the center point of cross-section BO that is perpendicular to the centerline. But it is not the 
center points of cross-section Bl and B2. In some cases non-perpendicular cross-sections result in a 
twisted or crooked centerline by changing the connecting order of center points. This correlation between 
orientation and center of a cross-section has been ignored in [26], which is also one of the main draw- 
backs in vessel tracking. The centerline also needs to be refined after being extracted. Refinement is an 
optimization process to approximate the centerline to the central axis, called the optimal centerline and 
also known as the good skeletonization [19]. Frangi et al. [1 1] introduce a refinement method in which a 
central axis and a vessel wall are modeled by B-spline surface and curve. The model is defined as snake 
and optimized by the energy field, which is directly defined by raw data (MRA). In this paper we present 
a general method to refine the centerline using a distance map. How to optimize the initial centerline in 
terms of its convergence is the main challenge. 

We define the centerline to be the center of vessel's walls (including lumen and plaque) rather than only 
its lumen. A separate paper will discuss a new technique of vessel segmentation in CTA. In this paper we 
use the distance map, named distance to boundary (DTB) volume, converted from binary data [24] to ex- 
tract and refine centerlines. 

In considering an optimization process, the rest of the paper is organized to answer the following ques- 
tions: 




• Existence: Does the optimal centerline exist in terms of the optimization goal? 

• Convergence: Does the algorithm converge to the optimal centerline if it exists? 

• Efficiency: How fast the algorithm will converge to the optimal centerline? 

In the next section we introduce the centerline definition for optimization, the principle of our optimization 
process, and the optimization iteration. The details of the optimization iteration based on the DTB volume 
are presented in section 3, including the computation of the reference frame, the center point, and the con- 
vergent process. Experiments and results are shown in section 4. We conclude our method in the last sec- 
tion. 

2 Motivation and Principle 

First of all, we need an unambiguous definition of a centerline in 3D for the purpose of optimization. 
Suppose that the vessel is a narrow tubular structure, which in general is a cylinder, the centerline can be 
regarded as the central curve axis of the cylinder. At each point of the central axis there is a cross-section 
that is perpendicular to the axis, i.e., the center of the cross-section is the center point of the centerline; 
the normal of the cross-section is the tangent of the centerline at this point. (In this paper when we men- 
tion cross-section, we refer to the cross-section that is perpendicular to the central curve axis.) Hence we 
define that the centerline consists of the centers of the cross-sections. We may define "center" in differ- 
ent ways, like geometrical center or physical centroid of the cross-section. 

The fundamental concept is the following recursive definition of a centerline (see Figure 3): 

• A centerline is a closure set of centers of the cross-sections of the object. 

• A cross-section is a cut plane that is perpendicular to the centerline. 

The fact is that we need a cross-section to compute a center point. But the position and orientation of a 
cross-section is defined by a segment of centerline, which is approximated or interpolated by a set of cen- 
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ter points. It is a "chicken & egg" problem. Fortunately we can obtain the initial centerline that we need 
to start the refinement process. 

We start with an initial centerline (which may be inaccurate); compute the cross-sections of the initial 
centerline; evaluate the center on each cross-section; then update the centerline by the center points 
evaluated. The new cross-sections will be computed by the updated centerline. This refinement process 
will not finish until the changes between successive loops is less than the accuracy required, i.e. when it 
converges to the optimal centerline. The iteration of this refinement process is shown in Figure 4. 

To adopt this definition, we still have the following questions not answered. 

• Is the iteration convergent? 

• If yes, does it converge to the optimal centerline, i.e. central axis? 

Supposing that a vessel segment is a cylinder, the cross-section at a center point is defined by the position 
(P) and the orientation (or tangent vector) (a) at this point. Thus, the area (5) of cross-sections within this 
segment is a function of P and a, i.e. S (P, a). The cross-section that is perpendicular to the centerline 
has the minimum area, i.e. min a { S (P, a) }. 

Assuming that the boundaries of a generalized cylinder are varied linearly within a small range, for sim- 
plicity we locate one boundary on the x-axis and another boundary on another line as shown in Figure 5. 
If these two boundaries are parallel, i.e., a uniform cylinder, it is obviously that the minimum area cross- 
section is perpendicular to the centerline, which locates at the middle of these two boundaries and is par- 
allel to the boundaries too. Considering that two boundaries are not parallel, the centerline is actually the 
angular bisector in terms of Blum's MAT. 

Supposing that P is a point on the angular bisector; the angle between an arbitrary oblique cross-section 
and the perpendicular cross-section is (3; the distance from \P to the boundaries is r; thus, the length of the 
oblique cross-section is 
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cos(/? + or) cos(/?-a) 
It is easy to prove that the shortest intersection line is found when /?= 0. 

Similar results are expected in the 3D case, see Figure 6. "When quadrilateral P u Q u QdP d rotates along X 
axis, SO (yellow cross-section) that is perpendicular to central axis (X axis) always owns the shortest in- 
tersection line compared to other cross-sections (blue cross-section) that are not perpendicular to the cen- 
tral axis. The area of cross-section is the integral of the area of all of the fans along the contours. Thus, 
the shortest intersection lines results in the minimum cross-sectional area. 

It shows us another fact of the centerline; each tangent vector of a centerline indicates the cross-section 
that has minimal cross-sectional area. The local minimum area ensures a unique convergent position. 
Therefore, the centerline refinement is an optimization process to find the orientation of minimum cross- 
sectional area within each segment, i.e. a cylinder with the centerline having n segments; S t is the cross- 
sectional area at segment i. 

Vimin(S,),i =1..» 

This concept of minimal cross-sectional area is reasonable in clinic practice. Imagine that we locate an 
oblique cut plane within a small segment of a vessel. There are many possible orientations and positions. 
The plane that we are most interested in is the one with minimum cross-sectional area in terms of steno- 
sis detection. This is the key idea of the optimization: to find the minimum cross-sectional area in each 
centerline segment. 

The pseudo code of one refinement process is shown below. The goal of refinement is to approximate the 
central axis by adjusting the points towards the cross-section centers, i.e. the optimal centerline (see 
Figure 7). 

In- Centerline: the centerline before refinement iteration 



Out-Centerline: the centerline after refinement iteration 

REFINEMENT STEPS: 

FOR (each Center-Point, on In-Centerline) 
o Calculate the Up and Normal (tangent) vectors of the cross-section at Center-Point, 
o Construct the Cross-Plane, in terms of Up, Normal vectors and the point Center-Point, 
o Fill the Cross-Plane, by scanning this plane in DTB volume 

o Detect the cross-section (we assume it is a general ellipse) on the Cross-Plane, and calculate the 

shape parameters of the ellipse, including area, long axis, short axis, etc. 
o Construct a local general cylinder with Cross-Plane, and its neighbors 
o Find the New Center-Point, and the New Tangentj 

o Put the New Center-Point, and the New Tangent into the Local Out-Centerline 
End FOR 

Gobal Approx. Centerline = Local Opt. Centerline.Least Square ApproximationQ ; 
Out-Centerline = Global Approx. Centerline. Discrete ResampledQ ; 
End Of Algorithm 

At each center point a local general cylinder is set up with the ellipse parameters extracted from its 
neighbors, see the red dashed line in Figure 7. The updated center point is limited within this local cylin- 
der. The new centerline is approximated globally by a least square cubic curve and re-sampled to a set of 
center points after one loop. 
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3 Convergence to Centerline 

3.1 Initialization 

Based on a DTB volume, we can extract an initial centerline using any centerline extraction algorithms or 
even via hand-drawing a piecewise linear centerline. Different centerline algorithms do not significantly 
affect the refinement, but might affect the computation time. The initial centerline may not be accurate 
but should be located within the object. (Actually, it is not strictly limited, as we will discuss in section 
3.4). In this paper, we use the CEASAR method [3] to create the initial centerline. 

3.2 Arclength Parameterization 

Optimization divides a centerline into a number of line segments, in which minimum cross-sectional area 
is evaluated. This division is done via parameterization of the initial centerline. The initial discrete cen- 
terline is first approximated by a cubic spline; we used a NURBS curve in this paper. Then, the approxi- 
mated curve is re-sampled equidistantly with a pre-defined arc-length {X - we use 2mm in this paper) to 
create a new discrete set of center points. Each re-sampled center point represents a small centerline 
segment of length A. The tangent vector of the centerline is the initial orientation of the cross-section at 
that point. 

3.3 Reference Frame 

A generalized cylinder is constructed by cross-sections, which are properly centered on the central axis. 
Assuming that vessels are not severely twisted, we can reconstruct a vessel by extruding a reference 
frame among cross-sections along its centerline. 

The reference frame consists of a reference point (P 0 - the position of the frame on the centerline) and 
three orthogonal axes (T 0 , B 0 , N 0 ) that define the orientation (see Figure 8). Tis the unit tangent vector of 
the centerline; B is the bi-normal vector and N is the principal normal vector. 




The Frenet frame is a convenient method to construct the reference frame. Unfortunately, it is undefined 
wherever the curvature vanishes, such as at points of inflection or along the straight sections of a curve. 
Also, the curvature vector can reverse direction on either side of an inflection point, inducing a violent 
twist in a progression of Frenet frames. This problem was discussed in several papers [4]. 

The alignment solution is a mechanism to minimize the torsion by using a small, local rotation along the 
curve. The first frame can be computed using curvature, as we do in the Frenet frame (if the local curva- 
ture vanishes, we select vector N one of vectors perpendicular to the tangent vector). Given the initial 
frame, the subsequent frames are computed by minimizing the torsion among its neighbors, as shown in 
Figure 8. First, a rotation axis and a rotation matrix is computed by using T 0 and T,. Then the old frame 
{P 0 , T 0 ) is rotated such that the T 0 aligns itself with the T,. This rotation creates a new N and B. By mov- 
ing the rotated frame to P,, we create a new frame (P,, T,) with the minimum torsion to P 0 . The details of 
rotation matrix and frame definition can be found in [4]. Because vessels are asymmetric, especially at 
the location of plaques, cross-section alignment with minimizing torsion is very important to ensure a 
correct local generalized cylinder. 

Each reference frame corresponds to a cross-section of a centerline. By aligning the x-axis of cross- 
section to the ./V vector and the y-axis of cross-section to the B vector, the cross-section is related to an 
oblique cut plane in space. This plane is filled in to the DTB volume; see Figure 9. Thus, a DTB cross- 
section is created as we mentioned in Section 2. 

3.4 Center Point 

According to the definition of a generalized cylinder, the center of a cross-section is the center point of 
the central curve axis, i.e. the optimal centerline. The center of a cross-section can be the geometric cen- 
ter or the physical centroid. A straightforward method to compute the center point is to find all of the 
boundary pixels in the cross-section, i.e. the contour, and calculate the center point by using the contour 




detected [26]. Here we introduce a more efficient way, a central moment method to estimate the center of 
a DTB cross-section. 

Suppose that a DTB cross-section is a 2D discrete function j{x, y), the y'th moment about zero is defined 
below, 



We can compute the x component /4 of the mean and y component /jy of the mean as 

(//,, jUy) = (m 0 , ,m, 0 ) 

{Mx,fJy) is the centroid point C in Figure 10, where point P is the center of the reference frame - reference 
point. 

Furthermore, we can define the central moments fly as below and use these moments to compute the use- 
ful descriptors of the local shape of a cross-section. 



where moments ju 2 o and JU02 are the variance of x and v, //// is the covariance between x and v. 

By finding the eigenvalues and eigenvectors of the covariance matrix, we can estimate the shape of a 
cross-section, including the short axis, the long axis, the eccentricity, the elongation, and the orientation 
of the shape, supposing it is an ellipse in general. 
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The covariance matrix is 
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In addition, we may notice that center point C does not directly relate to the reference point P. It means 
that the initial point may be located outside the vessel contour if the cross-section contains the vessel be- 
ing refined. 

3.5 Local Cylinder 

A local cylinder is constructed with the shape parameters on the current cross-section and its neighbors. 
In the ideal situation, the position P and the tangent vector T of the local central curve axis (see the 
dashed line in Figure 11) can be directly calculated if all the cross-sections are symmetric. In practice, we 
need an optimization process to approximate the central axis via the local minimal cross-sectional area. 

The optimization model that we use to search the minimum cross-sectional area is a spring model with a 
stochastic perturbation. Each center point is connected to its adjacent point by a spring. Instead of the 
displacement between two points, the difference between two orientations is the factor of the spring 
force, for example,/= k (1.0 - T 0 • T,). In addition each point is limited within its local cylinder, i.e. the 
equidistant re-parameterization arc-length (A) mentioned in Section 3.2. 

In Figure 11, cross-section i accepts spring forces from both cross-section i+ and i-, and vice versa. The 
stable orientation is defined by a weight summation of T it T-,+ and T { ., where the weight is the spring coef- 
ficient. In each iteration step, the cross-sectional orientations are adjusted by spring forces, then a sto- 
chastic perturbation vector is used to search for the local minimum area. Accordingly, the centerline is 
refined with the goal of minimum cross-sectional area constrained to the spring forces. 

3.6 Main Iteration 

The main iteration of centerline optimization consists of the steps below. 

• Parameterize the centerline with a spline curve: to create the discrete set of center points sampled 
equidistantly on the centerline, i.e. set up the initial position and orientation of the reference frames. 




• Adjust the orientation of each reference frame with the spring force from its adjacent frames, i.e. to 
modify the orientation of the reference frames. 

• Give stochastic perturbations to the orientation of the current frame to find the local optimized frame 
that has minimum cross-sectional area. 

• The center of the local optimized frame is taken as the refined center point. 

• Examine the convergence, if satisfied, the algorithm is finished. Otherwise, update the centerline for 
the next iteration. 

The criteria of convergence are the displacement of the center point and the deviation of the orientation 
(normal vector) of the reference frame. Although we use the minimum cross-sectional area to optimize 
the local center point, the sum of all cross-sectional areas cannot be taken as the global property of the 
optimum due to the facts below. 

• The reference frame is equidistantly positioned on the centerline. During optimization, center points 
are adjusted and the curve length of the centerline varies. Thus the number and the position of the 
reference frame may vary at each iteration step. 

• Since the position of the frame varies at each iteration step and the local cross-sectional area of the 
object is inconsistent, the local minimum cross-sectional area has no consistency among different 
iterations. 

Both the displacement of the center points and the deviation of the tangent vector of a centerline are 
taken as the factors of convergence. If both are less than a pre-defined threshold after certain iteration 
steps, we may think that the centerline is convergent. Both the maximum and average of the displacement 
and deviation are considered (see equations below, in which k is the iteration step, DS is the displacement 
and DV is the deviation of tangent vector, C is the updated center point, P is the position of a reference 
frame, Tis the updated tangent direction and N is the normal of a reference frame). 
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4 Experiments and Results 

In order to evaluate the algorithm, we used two different kinds of data sets, phantom data sets and clini- 
cal data sets. We use phantom data sets to evaluate the expected properties of the algorithm as well as its 
accuracy. The clinical data sets are used to evaluate the algorithm in practice, mainly for its reproducibil- 
ity. 

The first phantom test evaluates the effectiveness of the refinement. It may be argued that smoothing and 
approximation is enough for the centerline refinement (as we mentioned in Section 1). For this purpose, 
we generated an initial centerline in a pipe phantom using the CEASAR method [3]. Referring to the 
pseudo code of refinement in Section 2, we created two centerlines; one ignores the refinement steps and 
the other uses the refinement steps. Both centerlines are approximated by NURBS curve followed by two 
times of smoothing. The results are shown in Figure 12. Comparing image (a), in which only approxima- 
tion and smoothing are applied, and image (b), in which the full refinement process is applied, we ob- 
serve that the centerline processed by refinement is much smoother than the centerline processed only by 
approximation and smoothing. 

In addition we observe that the quality of the centerline is improved, i.e., not only the centerline becomes 
smooth, but it becomes more "central" if we compare the upper part of the centerline, especially the red 
line segment of the centerline. In fact the non-centered red line in image (a) is not one of the special er- 
rors created by the CEASAR method. It is a common problem using Dijkstra's shortest path searching 
algorithm. The searching usually considers two factors, "central" and "forward", in terms of costs. Dis- 
tance to boundary keeps the point "central" while searching orientation keeps the point moves forward 
from start point to end point. Cost is a compromise of these two factors. At the bending part of the pipe, 
the "shorter path" is a straight line. Image (a) shows us the effects of compromise between these two fac- 



tors. The sum of the cost of the red segment in Image (a) is smaller than the one in Image (b) since the 
former is shorter than the latter, which is the effect of searching orientation. In addition, we marked the 
biggest cross-sectional area by a green circle and the smallest one by a red circle. Since the pipe's diame- 
ter decreases progressively from top to bottom, the refined centerline in Image (b) shows more reason- 
able results. In another words, the tangential vector of the refined centerline is more accurate than the un- 
refined one. This test shows us the effectiveness of the refinement. 

The second phantom test is designed to evaluate the convergence of the refinement in different tubular 
cases. We designed three pipe phantoms with different diameters, one uniform pipe, one with a protuber- 
ance to simulate an aneurysm (called the aneurysm pipe), and one with a concavity to simulate a stenosis 
(called the stenosis pipe). The voxel size of the data is 1mm 3 . The reconstructed surfaces are shown in 
Figure 13. In order to see the convergence for pipes with different diameters, we set different diameters 
for these three pipes, uniform pipe of 3mm diameter, aneurysm pipe of 3.92mm diameter with an aneu- 
rysm of maximum diameter 7.37mm, and stenosis pipe of 4.78mm diameter with a stenosis of minimum 
diameter 1.91mm. The convergence is measured by the displacement and deviation of the tangential di- 
rection. The test results are shown in Figure 14, in which both the maximum and average deviations are 
displayed. We observe that the displacement and the tangential vector share a similar convergent pattern. 
It indicates that both the position and orientation of the centerline are convergent after refinement. As a 
matter of fact, we may check one of them to detect the convergence of the centerline in practice. 

Another visual way to examine the convergence of the centerline is using the Curve Planar Reformation 
(CPR) image [13]. Unlike the oblique cut plane, the CPR plane is a curved plane that scans through the 
volume along the centerline. The center point is the middle point of the scan line that is perpendicular to 
the centerline and the viewing direction. Thus the centerline is stretched and forms the central line of the 
image (see green lines in Figure 15). The CPR image quality is related to the smoothness of the center- 
line position and tangential vector. Compared with the images rendered with the initial centerline and the 




refined centerline in the uniform pipe, we observe that the refinement greatly improves the continuity and 
the accuracy. In addition, the aneurysm and stenosis are clearly displayed in the CPR images. 

The third phantom test is to evaluate the convergence at bifurcations. As we know, it is difficult to find 
the correct object region at bifurcations and the cross-section with minimum cross-sectional area might 
not be the correct cross-section. The solution is to limit the displacement of the center point at each re- 
finement step based on its neighbors. Due to the spring forces from its neighbors, the abrupt displace- 
ment at bifurcations will be adjusted and the center point will be drawn back to the appropriate position 
and orientation. Some examples of the refinement at bifurcations are shown in Figure 16. The left image 
is the case of two-branch bifurcation, in which the centerline goes from trunk to branch. The middle im- 
age is the case of three-branch bifurcation, in which three trunks meet at the bifurcation. The right image 
is the case mixing the cases of left and right images. All images show the high quality centerline at bifur- 
cations after refinement. 

In Figure 17 we show some cross-sections near the upper bifurcation along the centerline in the right im- 
age of Figure 16. Three middle cross-sections (image b, c and d) are enlarged and shown on the right side 
in Figure 17. In all images, the red spot is the center of the frame and it is the position of the center point 
before a refining iteration. The white spot is the center of the object cross-section and it is the center 
point after a refining iteration. As we discussed in the Section 3.4, the center is calculated by the centroid 
on DTB cross-section. As seen in the magnification of images b and d, the centroid remains centralized 
in the center of the main trunk rather than the geometrical center of the boundary contour used in [26]. 
We may observe some partial volume effects near the boundary of the objects on the cross-sections. But 
these boundary pixels have a very low DTB value. It has minimal effect on the centroid. 

There are two kinds of abnormal cross-section near bifurcation, one is called a suspended cross-section; 
the other is called an invalid cross-section. 




In the magnification of image c we observe that the centroid (white spot) is far from the center of the 
frame (red spot). If the displacement between the frame center and the centroid is bigger than the local 
cylinder arc length (A in section 3.2) or the short radius, which is calculated by the eigenvalues and ei- 
genvectors of the covariance matrix in section 3.4, we will mark this cross-section as a suspended cross- 
section. For a suspended cross-section, the center point will be shifted only the minimum radius dis- 
placement toward the centroid. In addition, the suspended cross-section has a low influence on its 
neighbors through the spring force. 

Another abnormal case is called an invalid cross-section when the object is intersected to the edges of the 
cross-section. It happens when the cross-section is severely tilted or the cross-section is too small. It also 
commonly happens near the bifurcation. For an invalid cross-section, we keep the initial center point, i.e., 
no displacement is made at the invalid cross-section. 

Two clinical data sets are used in this paper to evaluate the stability of the algorithm in practice, one is a 
coronal artery data set scanned by EBCT (electron beam computed tomography), and the other is an ab- 
dominal aorta data set scanned by MSCT (multi-slice CT). The main purpose of the clinical evaluation is 
to evaluate the reproducibility of vessel quantification. Two tests are designed. One is the stability of the 
refinement itself; the other is the stability caused by the initial centerline. 

In other words, the stability of the refinement process is the convergence of the algorithm. In the EBCT 
data set, we measure the Left Anterior Artery (LAD) and show the diameters along the LAD detected af- 
ter 10, 20, 30, 40 and 50 steps of refinement (see Figure 18). Because the diameter is directly measured 
from the cross section that is calculated by the position and orientation at each center point, the stability 
of the diameter shows the convergence from another point of view. The result shows that the refinement 
is convergent toward the centerline. Due to the resampling after refinement, the center points at each loop 
may be slightly shifted. Though we observe variation at bifurcation, the variation is less than 0.2mm. 
Considering the voxel size of 0.35 16mm 3 , this variation is acceptable. 




Another test is designed to evaluate the stability of the algorithm caused by the different initial center- 
lines. The stability of the measurement requires that the results done by different users should be repro- 
ducible. In practice, user's clicks for the start and end points of a centerline are always different. Thus, 
the initial centerlines are different too. We selected an abdominal aorta data set since the aorta has a rela- 
tively large cross-sectional area, which would more likely generate different initial centerline than in 
smaller vessels. In this test the user defines the aorta by two mouse clicks both more proximal and more 
distal (one at slice #70 the other is at slice #190) and the centerline algorithm is performed 12 times. 
Figure 19 demonstrates one of the results, a refined centerline and the corresponding CPR image. The 
measured results are listed in the table below, including the length, minimum and maximum area, mini- 
mum and maximum diameters, and the maximum curvature. We calculated the mean value, the average 
deviation and its percentage concerning the mean value. We observe that the deviation is less than 3%. 





Mean 


Average Deviation (AV) 


AV / Mean (%) 


Length (mm) 


132.8 


0.5185 


0.39 


Min-Area (mm 1 ) 


109.1 


2.2962 


2.10 


Max-Area (mm 2 ) 


235.3 


6.0740 


2.58 


Min-Diameter (mm) 


11.6 


0.0518 


0.45 


Max-Diameter (mm) 


22.2 


0.6617 


2.99 


Max-Curvature ( Degree/ mm) 


0.315 


0.0074 


2.35 



Finally we use a carotid CTA data set to do another clinical test to measure the efficiency, as shown in 
Figure 20. The 3D segmented carotid and its centerline is presented in the 3D volume-rendering image 
(upper image). We notice that the heavily calcified wall is segmented as part of the vessel as we men- 
tioned in Section 1 . In the middle image we can see two contours are shown, one red outer contour is the 



vessel contour, and the green inner contour is the lumen contour. It is clear that the vessel center point is 
not the lumen center point. The difference between the vessel area and lumen area is displayed on the 
lower MPR image. 

As an efficiency test, we recorded the computation cost on an Intel Pentium 4 - 2.8GHz CPU machine. 
The vessel is about 120mm long. The total time cost of vessel quantification is about 7.6 seconds, in 
which segmentation takes about 3.8sec, converting to DTB volume takes about 0.5sec, finding centerline 
takes about 0.5s, centerline refinement takes 2.2sec. The remaining 0.6sec is the computation time for an 
environment update, copying and saving results. The exporting quantification set consists of global data 
and cross-sectional data. Global data include vessel length, volume, minimum/maximum diameter, mini- 
mum/maximum area. Cross-sectional data include minimum/maximum diameter, lumen area, vessel area, 
wall contour, lumen contour, and tortuosity. 

5 Conclusion and Further Work 

Vessel quantification needs a reproducible and optimal centerline. Centerline optimization plays an 
important role in attaining a reproducible centerline. We have proposed a new centerline definition based 
on the minimum cross-sectional area and the optimization model constrained by the spring force among 
cross-sections. The center points are moved toward the centroid of cross-section at each iteration, i.e., the 
local central axis. This leads to the optimal centerline. 

Although we have tested more than 60 CTA data sets by using this algorithm, including different parts of 
the body, we expect further validation of our algorithms in the near future. In addition, we are investigat- 
ing the general centerline definition based on the minimum cross-sectional area, which may help us to ex- 
tend this vessel algorithm to other applications. 




Figure 1 Maximal inscribed sphere in 3D creates medial surface 




Figure 2 Correlation between orientation and center point on a cross-section 



Figure 3 A generalized centerline defined by the center points on each cross-section along the centerline 
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Figure 4 The centerline refinement process 




Figure 5 The shortest intersection line is perpendicular to the angular bisector 
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Figure 6 The minimum cross-sectional area (S) is the cross-section perpen- 
dicular to the central axis of a cylinder 




Figure 7 The convergent centerline 
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Figure 8 The subsequent frames are computed by rotation from its preceding frames 




Figure 9 Cross-section of DTB field 
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Figure 10 Centroid of the cross-section 




Input Centerline 



Figure 11 Local cylinder 





(a) (b) 

Figure 12 (a) Centerline smoothed by NURB approximation and smoothing filter, (b) Centerline 
refined by the optimization process. The green circle indicates the frame with the biggest cross- 
sectional area. The red circle indicates the frame with the minimal cross-sectional area. 
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Figure 13 Reconstructed surfaces of three phantoms. From left to right: uniform pipe of 3mm diameter, aneu- 
rysm pipe of 3.92mm diameter with an aneurysm of maximum diameter 7.37mm, and stenosis pipe of 4.78mm 
diameter with a stenosis of minimum diameter 1.91mm 





Figure 14 Deviation of displacement (upper) and tangent (lower) of three phantoms. From left to right: uni- 
form pipe, aneurysm pipe and stenosis pipe. (♦ : maximum deviation, □ : average deviation) 
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Figure 15 Curved Planar Reformation (CPR) images of the phantoms. From left to right, initial centerline 
in uniform pipe, refined centerline in uniform pipe, in aneurysm pipe and in stenosis pipe. Green line is the 
centerline stretched on the curved plane. 




Figure 16 Centerline refinement at phantom's bifurcations 
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Figure 18 Left anterior descending (LAD) artery diameters graph. The diameter ranges from 0.74mm to 
3.04mm. The average diameter is 1.63mm. Voxel size is 0.3516mm 3 . 
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Figure 20 Upper: 3D volume rendering of segmented carotid. Middle: one cross-section 
on the centerline. Lower: Cross sectional MPR of the carotid centerline. 
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Vascular Patent Claims 



1. Automatic vessel enhancement 

a. A fast multi-level vesselness filter method by using MIP- Volume pyramid. 

b. Automatic pre-CTA filter to enhance the vesselness response in CTA data 
sets 

c. Specific to body part either specified by the user or automatically 
determined 

d. Every voxel in dataset can be pre-assigned (automatically) with a 
probability of being a vessel 

e. Vessel probability determined by 3D shape, voxel intensity, n th derivatives 
of the intensity, texture, and any combination of these 

f. Enhancement used to improve 

i. Visualization of vasculature 

ii. Separation from non- vasculature such as bones, soft tissue 

iii. Segmentation of vasculature and other objects similar to vessels or 
tubes (spinal column) 

iv. Segmentation of non-vasculature such as bones, soft tissue 

v. Automatic extraction of vessels using basic knowledge of anatomy 
based on approximate direction, curvature, connectivity to other 
vessels 

1 . Use of this technique to extract major coronary vessels 
(e.g., left anterior descending artery (LAD), left circumflex 
artery (LCX), and right coronary artery (RCA) ) 

2. Use of this technique to extract major neck vessels such as 
common carotid arteries, internal carotid artery, external 
carotid artery, and vertebral arteries. 

3. Use of this technique to extract arteries of the brain such as 
the Circle of Willis and basilar arteries. 

4. Use of this technique to extract abdominal arteries such as 
the abdominal and thoracic aorta, renal arteries, common 
iliac arteries 

5. Use of this technique to extract peripheral arteries such as 
the common femoral arteries 

2. Vessel Segmentation with front propagation 

a. Automatic vessel seed point centralized: when moving mouse on any 
images, the mouse position is always positioned to the nearest center point 
if available. 

b. Fast vessel central axis tracked in vesselness volume gives users an 
overview to interrogate the vessel that is going to be segmented. 

c. A unique speed function defined with vesselness for front propagation. 

d. Fast vessel segmentation method by time-tag narrow band front 
propagation 

e. Segmented vessel includes both lumen and calcium. 

3. Centerline optimization 

a. A progressive optimization algorithm to refine a centerline 



b. A new centerline definition for refinement in terms of the minimum cross- 
sectional area. 

c. Calculate the lumen and wall contours on each cross-section, as well as 
other geometric information about these two contours. 

4. Automatic selection of processing protocol using meta-data supplied as part of the 
scan procedure (scanner supplied DICOM data fields) 

a. Based on a set of simple user-customizable rules that determine body part 

b. Based on body part, automatic body-part-specific segmentation and vessel 
enhancement can proceed. 

i. Automatic segmentation of heart 

1 . Removal of ribs, but leaving the spine and sternum for 
reference 

ii. Automatic segmentation of lungs 

c. Visualization parameters can then be automatically selected which may 
include 

i. Selection of 3D viewpoints 

1 . Designed to match standard hospital procedures such as 
cardiac, aortic, or brain catheterization. 

2. Any user-customizable set of viewpoints 

ii. Selection of a set of aforementioned 3D viewpoints that are 
automatically captured and saved to digital or film media 

iii. Selection of contrast/brightness setting or set of settings (called 
"window/level" in the parlance of radiology) specific to the body 
part 

iv. Selection of 3D opacity transfer function (or set of transfer 
functions) specific to the body part 

5. Curved MPR (CPR) on biconvex slab with rotation 

a. Using modified MIP projection 

b. Using modified x-ray projection 

c. Using adjustable diameter slab MIP projection 

6. Luminal MPR (CPR) on biconvex slab with rotation 

7. Display in 3D of the double-oblique cross-sectional slice or slab location 
(sometimes called "slice shadow") for easy cross-reference 

8. Synchronization of views containing a specific annotation so that the specific 
annotation desired is then visible in all applicable views 

9. Selections of an arbitrary plane for a double-oblique slice or slab view by drawing 
a line on any other view to create the plane. The new plane becomes the 
extrusion of the drawn line extended into the image. 

10. Selection of an arbitrary curved plane (an arbitrary curve extruded in one 
direction about a central 3D curve) created by drawing a curved on any other view 
to create the plane. 

11. Adjustment of a double-oblique view (arbitrary plane) by tilting the plane about 
the center of the image in any arbitrary direction 

a. Adjustment (tilting) about a set of known axes (e.g., horizontal, vertical, 
diagonal, or about image perpendicular axis) 

b. Translation of the center of the view by 




12. Electronic cleansing of hard plaques from arteries so that hard plaque becomes 
invisible in 2D and/or 3D views so that the blood lumen can be clearly seen 
without obstruction. 

a. The user-specifiable selection of the intensity threshold for calcification 

b. Partial volume effects are removed by Gaussian filter 

13. Selection and storage of multiple blood vessels that can be rapidly reviewed by 
selecting them one after another with all views updating with information about 
that vessel 

14. Display of blood lumen information graphs along the selected vessel on curved 
MPR and luminal MPR views 

a. Display of lumen area 

b. Display of calcification area on top of lumen area 

c. Display of minimum diameter 

d. Display of any varying parameter in graphical form synchronized 
alongside the vessel data 

e. Calculation of relative proportions (such as percent stenosis) based on the 
selection of two points along the vessel 

15. Simplified segmentation by placing a seed point at a desirable voxel location, 
computing some similarity or desirability measure based on nearby (or global) 
information around the selected location, then allowing user to interactively grow 
parts of the dataset that are similar to the selected location and nearby 

a. Allowing selection of more and less of the desired part based on a slider 
concept using distance along some scale as a metric to determine how 
much to include. This allows the user to select as small or as large as 
desirable the amount to include with ease. 

b. The preview of the selection in all applicable views 

c. The ability to add the current selection to already existing selections 

d. The determination of desirability for intensity data can be in proportion to 
the absolute difference relative to the intensity at the seed point 

e. The determination of desirability can be in proportion to the vessel 
probability measure (see earlier claims on vessel probability) 

f. The determination of desirability can be in negative proportion to the 
vessel probability measure (helpful for selecting non-vessel structures) 

g. The determination of desirability can be in proportion to a texture 
similarity measurement (e.g., using vector quantization to texture 
chararcteristics) 

h. The determination of desirability can be in proportion to shape-based 
similarity measurements (e.g., using curvature or n th derivatives of the 
intensity) 

i. The determination of desirability can be in proportion to some linear or 
non-linear combination of the above characteristics 

16. The ability to provide an endoluminal flight along the centerline of a vessel object 
(just like virtual colonoscopy but for blood vessels). 

a. With the display of hard plaque (calcification) and soft plaque in varied 
colors for differentiation from the vessel wall. 



b. With the ability to move back and forth along the centerline by direct 
manipulation of some linear mechanism 

i. By clicking or dragging a mouse along an overview of the entire 
vessel (like an overhead map) 

ii. By scrolling the mouse wheel to scroll along the centerline 

c. To interactively tilt the viewpoint about without leaving the centerline of 
the vessel 

17. The ability to perform selection (either for curved path generation, or seed point 
selection, or vessel endpoint selection, or 2D/3D localization) by clicking on an 
image and the selection of the point along the 3D line defined by the click point 
extruded into the image determined by 

a. The first point of intersection with the 3D object in which the voxel 
opacity or accumulated opacity reaches a certain threshold 

b. The middle point of entrance and exit of the 3D object as determined by a 
voxel intensity threshold 

18. The ability to select from a single view an area based on a single seed point 
deposit and to automatically compute the perimeter of the object and other 
particulars such as minimum diameter, maximum diameter, etc. 

a. With the included area of the object determined by an automatically 
derived threshold range 

b. With the included area of the object determined by a similarity measure of 
intensity, texture, connectivity, and derivatives of the intensity 

c. Also, the act of selection creates a set of annotations that describe the key 
' characteristics of the area automatically and displays these to the user 

19. Remove ribs from chest images and extract the heart region 

a. Lung region segmentation: thresholding or self-adaptive vector 
quantization method 

b. Determination of the center of heart: applying maximum principal 
component analysis method to determine the Long and the Short axis of 
the region that enclosed between two lung regions. The intersection point 
of those 2 axis is the center. 

c. Extract heart region: applying the radiation ray method. Sent ray from the 
center along any direction. Label the first voxel that a ray hits either the 
lung region or the volume boundary. Remove all voxels that are located 
beyond the first-hit voxel along all rays. The left region is the rough region 
of heart without ribs around it since the rib is surround the outside of the 
lung region. 

d. Extract accurate heart region: The bones in the rough heart region are 
chest plate and spine. Those bones can be removed by constrained region 
growing. 




